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Abstract 

The propagation of cosmic rays in the Earth's atmosphere is simulated. Calculations of the omnidi- 
rectional differential flux of neutrons for different solar activity levels are presented. The solar activity 
effect on the production rate of cosmogenic radiocarbon by the nuclear-interacting and muon compo- 
nents of cosmic rays in polar ice is studied. It has been obtained that the ^^C production rate in ice by 
the cosmic ray nuclear-interacting component is lower or higher than the average value by 30% during 
periods of solar activity maxima or minima, respectively. Calculations of the altitudinal dependence 
of the radiocarbon production rate in ice by the cosmic ray components are illustrated. 

1 Introduction 

Cosmogenic radionuclides are produced in the Earth's atmosphere and in matter on the Earth's surface in 
nuclear reactions caused by primary and secondary cosmic rays. The data on the cosmogenic radionuclide 
concentration in natural archives, such as tree rings and soil and ice layers, bear information on changes 
in the primary cosmic ray intensity. 

Radioactive isotope ^'^C is accumulated in polar ice together with atmospheric air precipitation in ice 
at the stage of glacier ice formation from firn. The nuclide production in reactions caused by secondary 
cosmic rays (in situ production) is another mechanism by which radiocarbon is accumulated in ice. 
Spallation reactions of oxygen atoms by neutrons of the cosmic ray nuclear-interacting component mainly 
contribute to the radiocarbon production in the upper ice layers. At depths larger than the depth of 
penetration of the cosmic ray neutrons (several meters), ^^C is produced in reactions caused by cosmic 
ray muons. 

The ^^C production rate in ice by secondary cosmic rays depends on the cosmic ray intensity. Geo- 
magnetic field variations do not affect galactic cosmic ray fluxes in polar regions. Thus, data on the in 
situ ^^C concentration in polar ice samples can be used as chronicles of solar activity variations in the 
past [1]. The present work studies the effect of solar activity variations on the radiocarbon production 
rate in ice by different cosmic ray components. The aim of this work is to quantitatively estimate the 
effect of the considered factor. 

2 Numerical calculations 

2.1 Galactic Cosmic Ray Spectrum 

At energies higher than several hundreds of MeV, galactic cosmic rays make the main contribution to the 
particle flux in the near Earth space. Protons and alpha particles account for about 90 and 10% of the 
total number of galactic cosmic ray particles, respectively. The particle fraction of heavier nuclei does 
not exceed 1%. The flux of galactic cosmic rays has a high degree of isotropy p|. 
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For particle energies lower than 10 GeV/nucleon, the flux of galactic cosmic rays in the near Earth 
space depends substantially on the solar activity level. The expression for the differential energy spectrum 
of i particles with mass and charge numbers Ai and Zi has the following form in the force field model [2| : 

T(F^-T (Fif^) E{E + 2m,c^) 

ME) - J,..5(i? + <J>.)(^^^^)(^^^^^2m,c2) (1) 

where Ji,Lis{E) is the differential energy spectrum of the i particles of galactic cosmic rays in the 
local interstellar space; E is the particle kinetic energy per nucleon; rUpC^ is the proton rest energy; 
$i = {eZi/ Ai)(j); e is the elementary charge; and (p is the parameter of cosmic ray modulation in the 
heliosphere. 

The following expression from [3] is used as a Ji,Lis{E) approximation: 

1 (F)-r P^^^'"''" (2) 

where parameter E corresponds to the particle kinetic energy in GeV/nucleon, Cj is the normalization 
factor, Cp = 1.9 x 10'^ (m^ s sr GeV)~^ for protons, and Chc = 9.5 x 10^ (m^ s sr GeV/nucleon) ~^ for 
helium nuclei. 

Usoskin et al. (2005) [1] presented the results of a reconstruction of the solar modulation parameter 
based on neutron monitor data for 1951-2004. The monthly average values of modulation parameter 
(f) vary from 0.3 to 2 GV, and the average (p value for the considered time interval is 0.69 GV. In our 
calculations, we selected the minimal and maximal values of the modulation parameter equal to 0.3 and 
1.2 GV, respectively. A parameter value of 0.3 GV corresponds to the monthly average value of the 
modulation parameter at a minimum of cycle 22, and 1.2 GV corresponds to the average annual value of 
this parameter at the cycle maximum 

The form of the energy spectrum of galactic cosmic ray particles at low energies in the force field model 
depends on the charge-to-mass number ratio Z/A. For nuclei with charge numbers Z > 2, Z/A ~ 1/2. 
The contribution of nuclei with Z > 2 to cosmic ray particle fiuxes in the Earth's matter was determined 
by multiplying the data obtained for ''He nuclei of galactic cosmic rays by coefficient k, which is the 
ratio of the nucleon flux caused by the nuclei of elements with Z > 2 to such a flux caused by galactic 
cosmic ray ^He nuclei. Based on the parameters of the energy spectra presented in [5] , it was found that 
k = 1.44. 

We ignored the geomagnetic cutoff of the cosmic ray differential energy spectrum, since we consider 
the production of cosmogenic ^^C in polar ice at high geomagnetic latitudes. 



2.2 Atmospheric Model 

In the numerical model, we assumed that the Earth radius and the height of the atmosphere are 6371 
and 100 km, respectively. It was also assumed that the atmosphere includes 0.755, 0.232, and 0.013 mass 
fractions of nitrogen, oxygen, and argon, respectively [6]. The atmosphere was modeled in the form of 
concentric layers with a constant air density within a layer. It was assumed that the thickness of each 
atmospheric layer is 15 g/cm^. The total thickness of the atmosphere was taken to be equal to 1034 
g/cm^. The dependence of the air density on altitude corresponds to the average air density at latitudes 
higher than 60°N according to the COSPAR data. In the calculations, we assumed that the Earth's 
surface layer is composed of II2O ice with a density of 0.917 g/cm^. 



2.3 Physical Model 

The numerical code for modeling cascades of cosmic ray particles was written based on the GEANT4 9.4 
simulation toolkit [7]. The processes of production, propagation, and interaction of baryons (neutrons, 
protons, short-lived baryons, and their antiparticles) , leptons (electrons, positrons, and muons), mesons 
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(pions and kaons), light nuclei, and gamma rays were taken into account in the calculations. The standard 
set of processes was used in modeling electromagnetic interactions [8| . Particles of the electron-photon 
component of secondary cosmic rays with energies lower than 10 MeV were not considered. Photonuclear 
and electronuclear processes, as well as decay processes, were taken into account in modeling. Low- 
and high-energy models based on the parametrization of experimental data were used to describe the 
processes of inelastic hadron scattering by nuclei. The Bertini intranuclear cascade model was used to 
model the processes of inelastic nucleon and meson scattering by nuclei at energies lower than 6 GeV. A 
model of the intranuclear binary cascade was used to describe inelastic interactions of light nuclei. The 
processes of capture of negative mesons and those of neutron capture and nuclear fission were taken into 
account in the calculations. To describe the interaction between neutrons and nuclei, we used the G4NDL 
3.14 package of cross sections. This package includes data from the ENDF-B VI, JENDL, and other data 
libraries for neutron energies lower than 20 MeV and the JENDL/HE data for energies varying from 20 
MeV to 3 GeV. 

2.4 Calculation of Particle Fluxes 

A flux of galactic cosmic ray particles with an energy higher than E through an area of the upper 
atmospheric boundary is 



Foi = j dVL COS e j dE'Jo,{E') = tt j dE'Joi{E') 



2-K E E 

where Joi{E) is the directed differential flux of galactic cosmic ray particles of sort i in the near-Earth 
space and 9 is the angle between the particle momentum vector and the direction toward the nadir point. 
The angular distribution of the number of particles that fall into the atmosphere through the area unit 
element satisfies the relationship dFoi/d cos 9 ~ cos 9. 

The omnidirectional differential flux of i particles with energy £^ in a target at depth z was calculated 
using the formula: 

Fop ^ I , , Fo. 



^^''^ NopAE^\cos9iy "No^AE^ \ COS 9u\ 

where the p and a subscripts mean proton and alpha particle, respectively; Nq is the number of particles 
of primary cosmic rays for which modeling was performed; the summation over / means summation with 
respect to all i particles that cross the specified z level in a target and have energies in the interval 
[E — AE/2; E + AE/2]; 9ii is the zenith angle of the ith particle; and k = 1.44. Particles with | cos 6ii\ < 
0.001 for atmospheric layers and | cos 9ii\ < 0.01 for ice were ignored in sum ^. The 1/ cos^ coefficient 
under the sign of the sum is used to calculate the omnidirectional particle flux. 

In each case, the number of primary particles, for which numerical simulation was performed, varied 
from 5 to 10 millions for protons and from 2 to 3 millions for alpha particles. The particle energy was 
chosen according to differential energy spectrum ([T]) and ([2]). Energies from 100 MeV to 1000 GeV 
for protons and from 100 MeV/nucleon to 1000 GeV/nucleon for alpha particles were considered. The 
calculations were performed at the Saint Petersburg branch of Joint Supercomputer Center of the Russian 
Academy of Sciences. 
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Figure 1 : Our calculations of the omnidirectional differential neutron flux at sea level for the conditions of 
high geomagnetic latitudes and moderate solar activity (a stepped line). The data from |9j are presented 
by a dot-and-dash line; the measurements performed by [lO], by a solid line with dots; the calculations 
made by [11], by a dotted curve; the calculations performed by [12j, by a dashed line. 

3 Results of numerical simulations 

3.1 Omnidirectional Differential Flux of Neutrons at Sea Level for the Conditions 
of Moderate Solar Activity 

Figure [T] illustrates calculations of the energy dependence of the omnidirectional differential neutron flux 
at sea level and moderate solar activity. Modulation parameter was taken equal to 0.69 GV in the 
calculations. For comparison. Fig. [T] presents the energy dependences of the differential neutron flux 
according to the data of the studies [9l [IHl [TH [T2]. 

The energy dependence of the differential neutron flux published in [9] resulted from an approximation 
of the available differential flux measurements. The data correspond to the New York location. 

The differential flux of neutrons from |10| resulted from a normalization of the measurements in the 
Northern Hemisphere for high geomagnetic latitudes, moderate solar activity, and sea level according to 
|10j . The measurement error is caused by the error of cross sections of interactions between high-energy 
neutrons and the matter that were used to calculate the spectrometer response function. According to 
pro], the error of the calculated response function values is about 10-15% for neutron energies higher than 
150 MeV and is lower for lower neutron energies, producing the measurement error of the same order of 
magnitude. 

The omnidirectional differential neutron fluxes obtained in jTT] resulted from simulations of cosmic 
ray propagation in the Earth's atmosphere. In fTT], numerical calculations were performed using the 
GEANT and MCNP program packages. In spite of the fact that we also used the GEANT program 
package in our calculations, the results of our modeling differ substantially from the data obtained in 
|llj . The reason of such discrepancy may be in different cross-section data used in both calculations. 
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In [12], numerical calculations were performed using the PHITS program, JENDL/HE cross sections 
were used to model the interactions between nucleons and nuclei. The result from [12] presented in Fig. 
[T] corresponds to the calculations performed ignoring the effect of the Earth's surface. The results of our 
calculations are in good agreement with the data obtained in [12j for neutron energies higher than 10 
MeV. 

At neutron energies 10-300 MeV, the calculations are in good agreement (within 30%) with the 
measurements [lOj. At neutron energies higher than 400 MeV, the difference between the calculated and 
experimental differential fluxes is substantial and reaches 200%. A similar difference with the experiment 
is observed for the calculated differential neutron flux from |12j . One of the possible causes of this 
difference consists in that the total JENDL/HE cross-sections used in both simulations are overstated at 
these energies. The difference in the cross-sections of the total interaction in 15% can result in a twofold 
difference in particle fluxes on five to six free paths. 

The differential flux of neutrons of secondary cosmic rays near the Earth's surface depends on the 
landscape and surface matter chemistry. A relatively small amount of hydrogen atoms in soil matter 
can substantially affect the neutron flux value. Disturbances in the neutron differential flux caused by 
water in soil decrease with increasing neutron energy: when the neutron energy is 10 MeV, the neutron 
flux disturbance is no more than 30% |12j . This phenomenon can be one of the causes of the differences 
between calculated differential fluxes and experimental data flO] at low energies (see Fig. [l]). 

3.2 Omnidirectional Differential Flux of Neutrons during Periods of Solar Activity 
Minima and Maxima 

Figure 2 illustrates the calculated energy dependence of the omnidirectional differential neutron flux for 
an atmospheric thickness of 660 g/cm^, which corresponds to an altitude of about 3.5 km. The calculation 
results for cosmic ray modulation parameters of 0.3, 0.69, and 1.2 GV are presented. According to the 
calculations, the average ratio of the energy spectrum values for the periods of solar activity minima 
and maxima is about 1.6 for neutron energies reaching several hundred MeV and is smaller for higher 
energies. It is interesting that the ratio of the neutron monitor count rates during periods of solar 
activity minima and maxima is about 1.4 for high geomagnetic latitudes and lower atmospheric layers 
[13j . The difference between the calculated and neutron monitor data can be partially explained by the 
sensitivity of neutron monitors to the high-energy part of the neutron spectrum and to the cosmic ray 
muon component. Cascade processes can also result from interactions between atmospheric matter and 
galactic cosmic rays with higher energies than it follows from nuclear interaction models. As a result, 
calculations give a larger amplitude of variations in the secondary cosmic ray fluxes. 

4 Production of ^^C in ice by cosmic ray components 

4.1 Production of ^"^C in ice by nucleons of cosmic rays at high geomagnetic latitudes 

Reactions of oxygen nuclei spallation by neutrons make the main contribution (98%) to the production 
of radiocarbon in ice by nucleons of the nuclear-interacting component of cosmic rays [T^ . The results of 
the work |14| were used to calculate the average radiocarbon production rate in polar ice by the nuclear- 
interacting component of cosmic rays. The main part of radiocarbon (80-90%) is produced in reactions 
induced by neutrons with energies 30-300 MeV. Based on calculations of the neutron differential flux in 
the atmosphere presented in the previous section, we can anticipate that the ^^C production rate in ice 
will be lower or higher than the average value by 30% during periods of solar activity maxima or minima, 
respectively. 
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Figure 2: Calculations of the omnidirectional differential neutron flux for an atmospheric thickness of 
660 g/cm^, for high geomagnetic latitudes, and three values of the solar modulation parameter: 0.3 (the 
upper curve), 0.69 and 1.2 GV (the lower curve). 



4.2 Production of ^ C in ice by cosmic ray muons 

Experimental data on the muon flux in water [15], as well as the data on the reaction cross-sections and 
the calculation results from [HI [T7|, were used to calculate the rate of ^^C production by muons. 

The production of the considered cosmogenic nuclide can result from inelastic interactions between 
high-energy muons and matter nuclei. Muons with energies close to the average muon flux energy mostly 
contribute to the nuclide production [16]. In ice nuclide is mainly produced by high-energy muons at 
depths of about ten meters of water equivalent. The average muon energy in matter is about 10 GeV at 
such depths [16]. According to our calculations, the ratio of omnidirectional fluxes of high-energy muons 
with energies higher than 5 GeV for periods of solar activity minima and maxima is 1.07 in the lower 
atmospheric layers. 

At the end of its motion, a negative muon is trapped onto the outer shell by one ambient atom and 
descends to a lower atom energy level. Then, the muon either decays into an electron and neutrino or 
reacts with a nucleus. In the latter case, the nuclear reaction can result in the production of a cosmogenic 
radionuclide. We consider that muons, the energy of which is lower than the energy of the differential 
muon flux maximum, i.e., about 0.5 GeV according to [18], are stopping muons. According to our 
calculations, the ratio of omnidirectional fluxes of low-energy negative muons with energies lower than 
0.3 GeV during periods of solar activity minima and maxima is 1.2-1.3 for atmospheric altitudes not 
higher than 4 km. 

Fluxes of cosmic ray muons in the atmosphere are less sensitive to a change in the solar activity level 
than those of nuclear-interacting component particles. 
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4.3 Production of ^^C in Ice at Different Altitudes above Sea Level 

This subsection presents the calculation results of the integral radiocarbon production rate in ice by cosmic 
ray components Qi depending on the glacier height; Qi = Pi{z)dz, where Pi{z) is the radiocarbon 
production rate in ice by cosmic ray component i at depth z. The value of the integral radiocarbon 
production rate determines the nuclide concentration in regions where polar ice is accumulated jl9j . 

The cosmogenic radionuclide production rate in ice is proportional to the omnidirectional integral 
flux of particles near the Earth's surface. According to the results presented in [201 EI], the nuclide 
production rate by neutrons of cosmic rays in ice at height h above sea level may be presented as 

where Q^q is the nuclide integral production rate in ice at height h above sea level; Q^q is such a rate at 
sea level; x is the thickness of the atmosphere in g/cuP above a specified level at height h, and Ajf™ is 
the attenuation length of the high-energy neutron integral omnidirectional flux at atmospheric altitudes 
varying from to h. According to our modeling, the average attenuation length of the high-energy 
neutron omnidirectional flux (neutron energies higher than 10 MeV) is A^*™ = 130 it 0.3, 134 it 0.3, and 
138 ± 0.4 g/cm^ at the level of one standard error for minimal, moderate, and maximal solar activity, 
respectively, and for atmospheric altitudes no higher than 4 km. According to the results presented in 
[2T] for high geomagnetic latitudes and altitudes reaching several kilometers, Aj^*™ = 130 -1-7-3 g/cm^, 
and the variation in the parameter during the 11-year solar cycle is about 4% [52], which confirms the 
results achieved by us in the present work. 

The muon flux attenuation length in the atmosphere is a function of the muon energy. A linear 
dependence is observed between the attenuation length of a vertical differential muon flux and particle 
momentum p3] • According to [21] , the attenuation length of an omnidirectional flux of negative muons 
with an energy lower than 0.3 GeV is 230-240 g/cm^ for high geomagnetic latitudes and altitudes reaching 
several kilometers. According to our calculations, the attenuation length of an omnidirectional flux of 
low-energy muons is about 250 g/cm'^ for atmospheric altitudes not higher than 4 km. The nuclide 
production rate in reactions of capture of negative muons in ice at height h above sea level is 



Qt- = Q^-exp 



/ 1034 - X 



V 



where Q^-, Q^- and x have the same meaning as in expression Q and A^*™ is the attenuation length 



= 240 g/cm^. 

The nuclide production rate by high-energy muons in ice at height h above sea level is 



of the omnidirectional flux of stopping muons; in these calculations, A|^_ = 240 g/cm 



where A^^" is the attenuation length of the omnidirectional flux of high-energy muons in the atmosphere. 
We assume that the value is equal to the attenuation length of a muon flux with energies of 5- 

10 GeV in the atmosphere. According to the results presented in ^3j and our calculations, A^j" ~ 
1000 g/cm^. Figure 3 presents the dependences of the rate of radiocarbon production by different cosmic 
ray components in ice on the glacier height at high geomagnetic latitudes and moderate solar activity. 
The total integral rate of radiocarbon production is also presented here. 

5 Conclusions 

We presented here the simulation results of the propagation of cosmic rays in the Earth's atmosphere for 
solar activity minima and maxima. We considered data on solar activity in the current epoch. According 
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Figure 3: Integral rates of C production in ice by cosmic ray components in reactions caused by cosmic 
ray neutrons Qn (solid thin curve) and high-energy muons Q^f (dot-and-dash curve) and in reactions of 
trapping negative muons Q^- (dashed curve); total production rate is represented by a thick solid 
curve. 

to the calculations, the ratio of fluxes of high-energy neutrons at minimal and maximal solar activity can 
reach 1.6 for the near-Earth atmospheric layers. At the same time, this ratio is less than 1.1 for high- 
energy muons. The contribution of the cosmic ray neutrons to the total integral radiocarbon production 
rate increases with increasing glacier height and is predominant for heights of more than 2 km. The ^'^C 
production rate in ice at high altitudes and high geomagnetic latitudes depends substantially on the solar 
activity level. 
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